xarray-spatial

User Guide

Generate terrain Hillshade Slope Curvature Aspect Bump map NDVI Mean Proximity Viewshed Zonal Statistics

Let's use datashader to render our images...

In [1]:
import numpy as np
import datashader as ds
from datashader.transfer_functions import shade
from datashader.transfer_functions import stack
from datashader.transfer_functions import dynspread
from datashader.transfer_functions import set_background
from datashader.colors import Elevation

Generate Terrain Data

The rest of the geo-related functions focus on raster data (or rasterized data, after a previous Datashader step that returns an Xarray object). To demonstrate using these raster-based functions, let's generate some fake terrain as an elevation raster:

In [2]:
from xrspatial import generate_terrain

W = 800
H = 600

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))
terrain = generate_terrain(canvas=cvs)

shade(terrain, cmap=['black', 'white'], how='linear')
Out[2]:

The grayscale value above shows the elevation linearly in intensity (with the large black areas indicating low elevation), but it will look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops:

In [3]:
shade(terrain, cmap=Elevation, how='linear')
Out[3]:

Hillshade

Hillshade is a technique used to visualize terrain as shaded relief, illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which is based on slope and aspect.

In [4]:
from xrspatial import hillshade

illuminated = hillshade(terrain)

shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')
Out[4]:

You can combine hillshading with elevation colormapping to convey differences in terrain with elevation:

In [5]:
stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),
      shade(terrain,     cmap=Elevation,         alpha=128, how='linear'))
Out[5]:

Slope

Slope is the inclination of a surface. In geography, slope is amount of change in elevation of a terrain regarding its surroundings.

Datashader's slope function returns slope in degrees. Below we highlight areas at risk for avalanche by looking at slopes around 38 degrees.

In [6]:
from xrspatial import slope

risky = slope(terrain)
risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)

stack(shade(terrain,     cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(risky,       cmap='red',              how='linear', alpha=200))
Out[6]:

Curvature

Curvature is the second derivative of the surface, or the slope-of-the-slope.

  • A positive curvature indicates the surface is upwardly convex at that cell.
  • A negative curvature indicates the surface is upwardly concave at that cell.
  • A value of 0 indicates the surface is flat.

Units of the curvature output raster are one hundredth (1/100) of a z-unit. The reasonably expected values of all three output rasters for a hilly area (moderate relief) can vary from -0.5 to 0.5; while for steep, rugged mountains (extreme relief), the values can vary between -4 and 4. Note that it is possible to exceed this range for certain raster surfaces.

In [7]:
from xrspatial import curvature

curv = curvature(terrain)

curv.data = np.where(np.logical_and(curv.data > 1e-06, curv.data < 2e-06), 1, np.nan)

stack(shade(terrain,     cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(curv,        cmap='green',              how='linear', alpha=200))
Out[7]:

Aspect

Aspect) is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.

Below, we look to find slopes that face close to North.

In [8]:
from xrspatial import aspect

north_faces = aspect(terrain)
north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,
                                          north_faces.data < 10), 1, np.nan)

stack(shade(terrain,     cmap=['black', 'white'], how='linear'),
      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),
      shade(north_faces, cmap=['aqua'],           how='linear', alpha=100))
Out[8]:

NDVI

The Normalized Difference Vegetation Index (NDVI) is a metric designed to detect regions with vegetation by measuring the difference between near-infrared (NIR) light (which vegetation reflects) and red light (which vegetation absorbs).

The NDVI ranges over [-1,+1], where -1 means more "Red" radiation while +1 means more "NIR" radiation. NDVI values close to +1.0 suggest areas dense with active green foliage, while strongly negative values suggest cloud cover or snow, and values near zero suggest open water, urban areas, or bare soil.

For our synthetic example here, we don't have access to NIR measurements, but we can approximate the results for demonstration purposes by using the green and blue channels of a colormapped image, as those represent a difference in wavelength similar to NIR vs. Red.

In [9]:
import xarray as xr

from xrspatial import ndvi

rgba = stack(shade(terrain, cmap=Elevation, how='linear')).to_pil()
r,g,b,a = [xr.DataArray(np.flipud(np.asarray(rgba.getchannel(c))))/255.0 
           for c in ['R','G','B','A']]

ndvi = ndvi(nir_agg=g, red_agg=b)
shade(ndvi, cmap=['purple','black','green'], how='linear')
Out[9]:

Bump

Bump mapping is a cartographic technique that can be used to create the appearance of trees or other land features, which is useful when synthesizing human-interpretable images from source data like land use classifications.

dsgeo.bump will produce a bump aggregate for adding detail to the terrain.

In this example, we will pretend the bumps are trees, and shade them with green. We'll also use the elevation data to modulate whether there are trees and if so how tall they are.

  • First, we'll define a custom height function to return tree heights suitable for the given elevation range
  • dsgeo.bump accepts a function with only a single argument (locations), so we will use functools.partial to provide values for the other arguments.
  • Bump mapping isn't normally a performance bottleneck, but if you want, you can speed it up by using Numba on your custom height function (from datashader.utils import ngjit, then put @ngjit above def heights(...)).
In [10]:
from functools import partial

from xrspatial import bump

def heights(locations, src, src_range, height=20):
    num_bumps = locations.shape[0]
    out = np.zeros(num_bumps, dtype=np.uint16)
    for r in range(0, num_bumps):
        loc = locations[r]
        x = loc[0]
        y = loc[1]
        val = src[y, x]
        if val >= src_range[0] and val < src_range[1]:
            out[r] = height
    return out

T = 300000 # Number of trees to add per call
src = terrain.data
%time trees  = bump(W, H, count=T,    height_func=partial(heights, src=src, src_range=(1000, 1300), height=5))
trees       += bump(W, H, count=T//2, height_func=partial(heights, src=src, src_range=(1300, 1700), height=20))
trees       += bump(W, H, count=T//3, height_func=partial(heights, src=src, src_range=(1700, 2000), height=5))

tree_colorize = trees.copy()
tree_colorize.data[tree_colorize.data == 0] = np.nan
hillshaded = hillshade(terrain + trees)

stack(shade(terrain,        cmap=['black', 'white'], how='linear'),
      shade(hillshaded,      cmap=['black', 'white'], how='linear', alpha=128),
      shade(tree_colorize,  cmap='limegreen',        how='linear'))

Wall time: 2.2 s
Out[10]:

Mean

The datashader.mean function will smooth a given aggregate by using a 3x3 mean convolution filter. Optional parameters include passes, which is used to run the mean filter multiple times, and also excludes which are values that will not be modified by the mean filter.

We can use mean to add a coastal vignette to give out terrain scene a bit more character. Notice the water below now has a nice coastal gradient which adds some realism to our scene.

In [11]:
from xrspatial import mean

LAND_CONSTANT = 50.0

water = terrain.copy()
water.data = np.where(water.data > 0, LAND_CONSTANT, 0)
water = mean(water, passes=50, excludes=[LAND_CONSTANT])
water.data[water.data == LAND_CONSTANT] = np.nan

stack(shade(terrain,    cmap=['black', 'white'], how='linear'),
      shade(water,      cmap=['aqua',  'white']))
Out[11]:

Full scene

We've now seen several of datashader's geo helper functions for working with elevation rasters.

Let's make a full archipelago scene by stacking terrain, water, hillshade, and tree_colorize together into one output image:

In [12]:
stack(shade(terrain,       cmap=Elevation,          how='linear'),
      shade(water,         cmap=['aqua','white']),
      shade(hillshade(terrain + trees),   cmap=['black', 'white'], how='linear', alpha=128),
      shade(tree_colorize, cmap='limegreen',        how='linear'))
Out[12]:

Proximity

The datashader.spatial.proximity function operates on a given aggregate to produce a new distance aggregate based on target values and a distance metric. The values in the new aggregate will be the distance (according to the given metric) between each array cell (pixel) and the nearest target value in the source aggregate.

A powerful feature of proximity is that you can target specific values in the aggregate for distance calculation, while others are ignored. Play with the target_values parameter below and see the difference of using target_values=[1,2,3] vs. target_values=[2] vs. target_values=[3]

Load data and create ds.Canvas
In [13]:
from xrspatial import proximity
import pandas as pd


df = pd.DataFrame({
   'x': [-13, -11, -5,4, 9, 11, 18, 6],
   'y': [-13, -5, 0, 10, 7, 2, 5, -5]
})

cvs = ds.Canvas(plot_width=W, plot_height=H,
                x_range=(-20, 20), y_range=(-20, 20))
Create Proximity Aggregate
  • Use Canvas.points to create an xarray.DataArray
  • Calculate proximity to nearest non-nan / non-zero elements using datashader.spatial.proximity
In [14]:
points_agg = cvs.points(df, x='x', y='y')
points_shaded = dynspread(shade(points_agg, cmap=['salmon',  'salmon']),
                          threshold=1,
                          max_px=5)
set_background(points_shaded, 'black')
Out[14]:
Create proximity grid for all non-zero values
In [15]:
proximity_agg = proximity(points_agg)

stack(shade(proximity_agg, cmap=['darkturquoise', 'black'], how='linear'),
      points_shaded)
Out[15]:
In [16]:
line_agg = cvs.line(df, x='x', y='y')
line_shaded = dynspread(shade(line_agg, cmap=['salmon',  'salmon']),
                          threshold=1,
                          max_px=2)
set_background(line_shaded, 'black')
Out[16]:
In [17]:
line_proximity = proximity(line_agg)
stack(shade(line_proximity, cmap=['darkturquoise', 'black'], how='linear'),
      line_shaded)
Out[17]:
Transform Proximity DataArray

Like the other Datashader spatial tools, the result of proximity is an xarray.DataArray with a large API of potential transformations.

Below is an example of using DataArray.where() to apply a minimum distance and maximum distance.

In [18]:
where_clause = (line_proximity > 1) & (line_proximity < 1.1)
proximity_shaded = shade(line_proximity.where(where_clause), cmap=['darkturquoise', 'darkturquoise'])
proximity_shaded = set_background(proximity_shaded, 'black')
stack(proximity_shaded, line_shaded)
Out[18]:

Viewshed

The datashader.spatial.viewshed function operates on a given aggregate to calculate the viewshed (the visible cells in the raster) for the given viewpoint (observer) location.

The visibility model is the following: Two cells are visible to each other if the line of sight that connects their centers does not intersect the terrain. If the line of sight does not pass through the cell center, elevation is determined using bilinear interpolation.

Simple Viewshed Example
  • The example below creates a datashader aggregate from a 2d normal distribution.
  • To calculate the viewshed, we need an observer location.
  • This location is indicated by the orange point in the upper-left of the plot.
In [19]:
from xrspatial import viewshed

OBSERVER_X = -12.5
OBSERVER_Y = 10

canvas = ds.Canvas(plot_width=W, plot_height=H,
                   x_range=(-20, 20), y_range=(-20, 20))

normal_df = pd.DataFrame({
   'x': np.random.normal(.5, 1, 10000000),
   'y': np.random.normal(.5, 1, 10000000)
})
normal_agg = canvas.points(normal_df, 'x', 'y')
normal_agg.values = normal_agg.values.astype("float64")
normal_shaded = shade(normal_agg)

observer_df = pd.DataFrame({'x': [OBSERVER_X], 'y': [OBSERVER_Y]})
observer_agg = canvas.points(observer_df, 'x', 'y')
observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),
                            threshold=1, max_px=4)

normal_illuminated = hillshade(normal_agg)
normal_illuminated_shaded = shade(normal_illuminated, cmap=['black', 'white'], 
                                  alpha=128, how='linear')

stack(normal_illuminated_shaded, observer_shaded)
Out[19]:
Calculate viewshed using the observer location
In [20]:
# Will take some time to run...
%time view = viewshed(normal_agg, x=OBSERVER_X, y=OBSERVER_Y)

view_shaded = shade(view, cmap=['white', 'red'], alpha=128, how='linear')

stack(normal_illuminated_shaded, observer_shaded, view_shaded)

Wall time: 7.07 s
Out[20]:
Viewshed on Terrain
  • Let's take the example above and apply it to our terrain aggregate.
  • Notice the use of the observer_elev argument, which is the height of the observer above the terrain.
In [21]:
from xrspatial import viewshed


cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))
terrain = generate_terrain(canvas=cvs)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')

illuminated = hillshade(terrain)

OBSERVER_X = 0.0
OBSERVER_Y = 0.0

observer_df = pd.DataFrame({'x': [OBSERVER_X],'y': [OBSERVER_Y]})
observer_agg = cvs.points(observer_df, 'x', 'y')
observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),
                            threshold=1, max_px=4)

stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),
      terrain_shaded,
      observer_shaded)
Out[21]:
In [22]:
%time view = viewshed(terrain, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=100)

view_shaded = shade(view, cmap='fuchsia', how='linear')
stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),
      terrain_shaded,
      view_shaded,
      observer_shaded)

Wall time: 2.26 s
Out[22]:

The fuchsia areas are those visible to an observer of the given height at the indicated orange location.

Zonal Statistics

Zonal statistics allows for calculating summary statistics for specific areas or zones within a datashader aggregate. Zones are defined by creating an integer aggregate where the cell values are zone_ids. The output of zonal statistics is a Pandas dataframe containing summary statistics for each zone based on a value raster.

Imagine the following scenario:

  • You are a hiker on a six-day-trip.
  • The path for each day is defined by a line segement.
  • You wish to calculate the max and min elevations for each hiking segment as a Pandas dataframe based on an elevation dataset.
In [23]:
from datashader.colors import Set1

cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))

terrain = generate_terrain(canvas=cvs)
terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')

illuminated = hillshade(terrain)
illuminated_shaded = shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')

zone_df = pd.DataFrame({
   'x': [-11, -5, 4, 12, 14, 18, 19],
   'y': [-5, 4, 10, 13, 13, 13, 10],
   'trail_segement_id': [11, 12, 13, 14, 15, 16, 17]
})

zones_agg = cvs.line(zone_df, 'x', 'y', ds.sum('trail_segement_id'))
zones_agg.values = np.nan_to_num(zones_agg.values, copy=False).astype(np.int)
zones_shaded = dynspread(shade(zones_agg, cmap=Set1), max_px=5)

stack(illuminated_shaded, terrain_shaded, zones_shaded)
Out[23]:
In [24]:
from xrspatial import zonal_stats

zonal_stats(zones_agg, terrain)
Out[24]:
mean max min std var
11 1321.919879 1428.384240 1228.065037 53.409131 2852.535316
12 1489.600166 2399.556289 1220.694206 310.177594 96210.140061
13 1711.632114 2430.430303 1368.922124 278.369664 77489.669636
14 1355.661087 1368.952389 1346.856316 6.249734 39.059170
15 1381.950485 1443.856734 1319.891060 38.918956 1514.685125
16 1486.314866 1709.222348 1367.673562 120.950684 14629.068060
Calculate custom stats for each zone
In [25]:
custom_stats = dict(elevation_change=lambda zone: zone.max() - zone.min(),
                    elevation_min=np.min,
                    elevation_max=np.max)

zonal_stats(zones_agg, terrain, custom_stats)
Out[25]:
elevation_change elevation_min elevation_max
11 200.319203 1228.065037 1428.384240
12 1178.862083 1220.694206 2399.556289
13 1061.508179 1368.922124 2430.430303
14 22.096073 1346.856316 1368.952389
15 123.965674 1319.891060 1443.856734
16 341.548786 1367.673562 1709.222348

Here the zones are defined by line segments, but they can be any spatial pattern, and in particular can be any region computable as a Datashader aggregate.

References